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UNITS OF MEASURE 


In a prepared statement presented on August 5, 1965, to the 
U. S. House of Representatives Science and Astronautics Committee 
(chaired by George P. Miller of California), the position of the 
National Aeronautics and Space Administration on Units of Measure 
was stated by Dr. Alfred J. Eggers, Deputy Associate Administrator, 
Office of Advanced Research and Technology: 

"In January of this year NASA directed that the international 
system of units should be considered the preferred system of units, 
and should be employed by the research centers as the primary 
system in all reports and publications of a technical nature, except 
where such use would reduce the usefulness of the report to the 
primary recipients. During the conversion period the use of cus- 
tomary units in parentheses following the SI units is permissible, 
but the parenthetical usage of conventional units will be discontinued 
as soon as it is judged that the normal users of the reports would 
not be particularly inconvenienced by the exclusive use of SI units." 

The International System of Units (SI Units) has been adopted 
by the U. S. National Bureau of Standards ( see NBS Technical News 
Bulletin, Vol. 48, No. 4, April 1964) . 

The International System of Units is defined in NASA SP-7012, 
"The International System of Units, Physical Constants, and 
Conversion Factors," which is available from the U. S. Government 
Printing Office, Washington, D. C. 20402. 

SI Units are used preferentially in this series of research re- 
ports in accordance with NASA policy and following the practice of 
the National Bureau of Standards. 
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PREFACE 


In 1955, the team which has become the Marshall 
Space Flight Center (MSFC) began to organize a re- 
search program within its various laboratories and 
offices. The purpose of the program was two-fold: 
first, to support existing development projects by re- 
search studies and second, to prepare future develop- 
ment projects by advancing the state of the art of 
rockets and space flight. Funding for this program 
came from the Army, Air Force, and Advanced Re- 
search Projects Agency. The effort during the first 
year was modest and involved relatively few tasks. 
The communication of results was, therefore, com- 
paratively easy. 

Today, more than ten years later, the two-fold 
purpose of MSFC's research program remains un- 
changed, although funding now comes from NASA Pro- 
gram Offices. The present yearly effort represents 
major amounts of money and hundreds of tasks. The 
greater portion of the money goes to industry and uni- 
versities for research contracts. However, a sub- 
stantial research effort is conducted in house at the 
Marshall Center by all of the laboratories. The com- 
munication of the results from this impressive re- 
search program has become a serious problem by 
virtue of its very voluminous technical and scientific 
content. 

The Research Projects Laboratory, which is the 
group responsible for management of the consolidated 
research program for the Center, initiated a plan to 
give better visibility to the achievements of research 
at Marshall in a form that would be more readily us- 
able by specialists, by systems engineers, and by 
NASA Program Offices for management purposes. 

This plan has taken the form of frequent Research 
Achievements Reviews, with each review covering one 
or two fields of research. These verbal reviews are 
documented in the Research Achievements Review 
Series. 


Ernst Stuhlinger 

Director, Research Projects Laboratory 

These papers presented January 6, 1966 
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RESEARCH IN COMPUTATIONAL MATHEMATICS AND LANGUAGES 


By 

C. L. Bradshaw’ 


SUMMARY 

QO3°0 

The problems associated with effectively relating 
machine languages and problem oriented language for 
efficient use of computers is discussed. The need 
for improvements in the man-machine relationship 
and economic improvements of trade-offs that can be 
achieved between presently developed machine lan- 
guages to broaden the participation and usefulness of 
the computer in space programs is emphasized. 

The impact of standardization on computer utili- 
zation, the consequent reduction in redundant effort 
and relief from the continual need for the reformula- 
tion of the problem is shown as an important objective 
of the Computation Laboratory and some achievements 
in this area are discussed. 


I. INTRODUCTION 


Computers have become an essential tool in the 
research and development programs of our nation. 
They have also become a very expensive and sensi- 
tive item in our nation’s budget. Research into 
computational mathematics and languages can lead to 
a more effective use of this most important tool. 

This research effort can attack the overall problems 
on five main fronts: 

1. Improving the mathematics involved in obtain- 
ing a computer solution to an engineering or scientific 
problem, 

2. Improving the computer programming lan- 
guages used in problem solution, 

3. Obtaining a more effective use of computing 
hardware and software as relates to specific classes 
of problems, 

4. Development of more efficient computer 
hardware, and 


5. Improving the man-machine relationship as 
relates to automatic computing devices. 

MSFC has done or sponsored considerable re- 
search in areas which have direct application to the 
areas mentioned above. This survey will show the 
MSFC efforts in these areas. This presentation will 
be followed by two papers which will be more explicit 
in two of these areas. 


II. COMPUTATIONAL LANGUAGES 


As mentioned earlier, computers are now being 
used to solve many diverse problems in the fields of 
engineering, science, and business. The fast com- 
puting speed and large internal memory of general 
purpose computers are valuable assets to those who 
prepare the problem solution. However, the language 
of computers is a sequence of numbers which usually 
is reduced by the machine to a sequence of ones and 
zeros, and this machine language is generally foreign 
to the problem solution prepared by an expert in a 
specialized area. Thus, there is a gap between the 
language of a problem and the language of a machine 
for solving it. The seriousness of this gap is in- 
tensified by the fact that there are almost as many 
machine languages as there are kinds of computers. 
Therefore, work performed in machine language at 
one site is often of little value elsewhere if the 
machines are different. 

To reduce this gap, many programs have been 
prepared so that on the one hand they can be under- 
stood by computers and on the other hand they will 
accept as input a higher level language which is closer 
to the problem and is called the source language . The 
techniques to implement programs which accept 
source languages fall into two overlapping categories. 
The simplest is to have the computer interpret state- 
ments of the source language and process the intent 
when it is recognized. The more popular technique 
is to translate the source program into an object pro- 
gram which is either in machine language or closer 
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to machine language. In the latter technique, the 
source program may have several equivalent repre- 
sentations since the object language from one trans- 
lation may be the source language for another. 

Pragmatically, a language is defined in terms of 
what a processor (whether it be called compiler, 
translator, or assembler) will recognize correctly 
for the user. In other words, a computer language 
is directly or indirectly the language for input to a 
computer. Consequently, standardization with a 
group of users is necessary to determine what is 
mutually agreed to as being the language. 

When the users of a language are allowed to 
participate in the formulating process, a more mature 
version of the language can be produced more quickly. 
The definition of ATOLL II (Automatic Test Oriented 
Launch Language) , which we will now describe, il- 
lustrates how the changing needs of the user can be 
incorporated into the language design during definition 
to avoid producing numerous languages. In such an 
environment, the trade-offs between features desired 
and the implementation cost can be evaluated realis- 
tically, 

ATOLL II is a problem -oriented language for 
real-time launch vehicle testing. The language is 
structured like FORTRAN (Formula Translation) and 
includes, in addition, real-time test-oriented state- 
ments, a more elaborate data description capability, 
a limited ability to manipulate symbolic or other non- 
numeric values, and an ability to include inline sym- 
bolic coding. 

The language provides the capability to manipu- 
late both the ground support equipment and the launch 
vehicle. It provides for real-time delays, for con- 
trol based on maintaining a sequence of events where 
event execution is time related with respect to pre- 
vious functions. 

The language is open ended in that the user may 
define what appears to the user to be additional source 
statement types. Thus, the compiler can be adapted 
to the problem area as new specifications are pre- 
pared in the language. 

The language is block structured to permit dy- 
namic allocation of variable and temporary storage. 
This feature, in combination with a provision to auto- 
matically segment a test into independent programs, 
assures that the object space required for test pro- 
grams can be kept very small. 

The language has been designed to satisfy many 
of the requirements and desires of management, the 


launch system engineer, and the computer program- 
mer. According to dynamic definition techniques 
which have been developed, this language has evolved 
rapidly over the past year. ATOLL II is fully docu- 
mented and available for study. 

We now would like to mention the problem of 
language translation. Since much programming ef- 
fort has already been expended in languages which 
are now obsolete, or for which processors are not 
readily available, there is interest in the capability 
to translate a program into another language which 
is available. Also, it is desirable to minimize the 
number of different languages a user must learn. 

Many super-processors have been proposed 
which are "machine independent." The purpose of 
such a processor is to allow preparation of compilers 
for classes of languages rather than for only a spe- 
cific one. Until recently, insufficient information 
was available to determine what classes are suffi- 
ciently defined for implementing in such general 
terms. However, some efforts have been quite suc- 
cessful in the areas of assemblers and context-free 
languages. How to proceed with context-sensitive 
languages is not yet clear, although preliminary 
efforts in this direction are being made. 

Some of the most dramatic developments in soft- 
ware have been seen in operating systems, or so- 
called control programs. Control programs are 
being given the tasks of handling computer interrupts, 
doing bookkeeping on jobs, servicing remote termi- 
nals on a priority basis, scheduling memory and 
computer time, editing and merging of programs, 
and total data management. As techniques of mech- 
anizing such tasks are developed, the user is freed 
from meticulous operations and his turn-around time 
is shortened so that his time can be used more pro- 
ductively. 

One problem area which has experienced many 
attempts but little success is that of preparing lan- 
guage processors. Early attempts such as Jovial 
and Neliac have proved to be educational but econom- 
ically unfavorable. Recent developments in assembly 
level languages have improved expressibility, strati- 
fied the control of symbol expressions, incorporated 
list structures, and refined recursive macro capa- 
bility with conditional parameter substitution. Yet, 
more work is needed to fill the gap between the kinds 
of languages which are easily implemented and the 
kinds of statements which users in specialized areas 
find most appropriate for the problem at hand. 

An example of a specialized area which has justi- 
fied the development of a new language is trajectory 
programming, which is discussed in the next section. 
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1 1 1. THE MSFC TRAJECTORY LANGUAGE 

In the past, MSFC scientific programmers have 
worked individually with the engineer in the develop- 
ment of programs which were designed, programmed, 
checked out, and documented specifically for that 
engineer’s need. This relationship of programmer 
to engineer has proved effective because of the nature 
of past problems and the limited language and sys- 
tems capabilities. Specifically, in the past, many 
programs were large and involved with long produc- 
tion lifetimes. More recently, however, desired 
programs cover a wider range of applications with 
more limited use. This fact makes it imperative 
that our scientific programmers produce and main- 
tain many more programs. Also, software advances 
have been made which open the door to a more general 
and sophisticated approach to the trajectory applica- 
tions area. There seems to be no choice but to re- 
evaluate our overall procedures and optimize where 
possible. As the result of careful investigation, it 
has been determined that the function of setting up 
and maintaining trajectory programs can and must be 
optimized to a maximum reasonable level. At this 
point, the maximum level of optimization or automa- 
tion cannot be ascertained; however, some optimiza- 
tion can be realized. This is possible when one 
realizes that the entire area of trajectory computa- 
tion, when taken as a whole, is a set of associated 
problems with many elements in common. 

The form that a problem may acquire in the pro- 
cess of being prepared for computation will vary 
widely within the range of imagination, experience, 
and other resources possessed by individuals who 
perform this task. For this reason, it is frequently 
difficult and time consuming for one person to use or 
become familiar with a program that was written by 
someone else. It is also a time consuming task for 
a programmer to modify his own program. Detailed 
documentation relieves this problem to some extent, 
but the more documentation there is to be studied, 
the more complicated the task becomes. It is neces- 
sary, or at least desirable, to restrict the general 
overall structure of programs to conform as nearly 
as possible to a general well defined standard model. 
As a result of the language and systems improve- 
ments and increased workloads, it was imperative 
that research be done to establish a faster, less ex- 
pensive, and more useful service by developing a 
trajectory oriented programming system. 

With the increasing demands for faster results 
from man and the computer, it has become obvious 
that these demands cannot be met with present re- 
sources. With the cost of manhours increasing and 


the cost of computing machine time decreasing, more 
burden must be placed on the computer by the use of 
more problem oriented systems. The trajectory pro- 
grammer must be enabled to do a better job in a 
shorter time and at less cost to meet these demands. 

We answer the question, "How can a trajectory 
oriented programming system help eliminate some 
of the effort required by the man in the man -to - 
machine cycle?" as follows: 

1. Construct the programs in modular form. 
Programs could consist of elements called modules 
and can be thought of as building blocks for many 
programs. This will eliminate repetitive efforts 
since modules may be interchanged to create com- 
pletely different programs without major reprogram- 
ming effort. 

2. Standardize nomenclature. To implement the 
modular concept, standardization is necessary. It 
enables the programmer and the user to communicate 
in well defined terms, eliminating confusion in defi- 
nition of coordinate systems, mathematical models 
and units. Documentation of work done will be more 
effective and meaningful. Programs written by other 
programmers will be easier to interpret and under- 
stand. 

3. Standardize organization of programs. The 
user and programmer can communicate at a common 
level. The programmer would be free to do more 
useful and creative work in other areas. New users 
and new programmers can become familiar with pro- 
grams and trajectory concepts earlier. Interchanging 
of programs will be easier. Programs will be easier 
to modify and maintain simply by changing and modi- 
fying only the necessary modules. Programs will be 
easier to evaluate since the programming effort will 
be isolated from system functions. Organization of 
the problem will be simplified since much of the 
logic will be handled by the system preprocessor. 

We next look at the impact on utilization. 

A problem oriented language and system will 
enable the programmer to drastically cut the time 
required in setting up and maintaining a trajectory 
related program. This savings results from the 
programmer being able to use precoded subroutines 
and sub-programs as the need arises. These pre- 
coded elements will be fully documented and com- 
pletely checked out beforehand, thus freeing the 
programmer from these routine tasks. 

A library will be established for the programmer 
and the user. Therefore, the engineer will be relieved 
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of the never ending problem of reformulation and re- 
checking of requests and proposals. As routines are 
developed and pooled by users, much redundant effort 
will be eliminated. 

Better documentation of work will be provided by 
a standardization of program logic. 

We now look at the difficulties in developing a 
trajectory system of this type. To be widely useful 
and accepted, a trajectory system must satisfy the 
requirements of all the users. Each user usually has 
a special interest which emphasizes certain areas 
more than others; for instance, the engineer whose 
primary concern is trajectory optimization puts a 
different emphasis on the guidance package than will 
the engineer whose primary concern is simulating 
the on-board guidance computer. In actual practice, 
it is difficult to satisfy the needs of the many users 
with one trajectory system. Even though difficult, 
sufficient research has been done to demonstrate the 
feasibility of further research of such effort. 

We must recognize the research that must be 
carried on in the development of a problem oriented 
system. 

1. Problem areas must be analyzed to determine 
the feasibility of developing and implementing a tra- 
jectory language, in our case, to see if the trajectory 
area is sufficiently large to justify such an effort. 

We have determined that it is large enough. 

2. Work done by other laboratories and installa- 
tions must be investigated to determine if they are 
doing work in this area. If so, we must establish 
why the effort , amount of effort , state of development , 
and evaluation of their effort relative to our needs. 

3. Our own needs, present and future, must be 
investigated considering customer contacts, past 
requirements, survey of existing programs, what 
future needs are expected to be. 

4. Methods must be developed for providing the 
following: 

a. Program logic - to allow the programmer 
to direct the flow of the program in extremely in- 
volved logical paths. 

b. Events - all complex trajectory programs 
involve some type of events or interrupts; much plan- 
ning and analysis are required prior to writing a 
program. Examples of events are engine cutoffs, 
weight drops, high-q and tilt arrests. 


c. Integration - much research needs to be 
done in this area to allow the user to select integra- 
tion schemes to give the desired accuracy for his 
problems. 

d. Input and Output - the present laborious 
effort of writing long lists and involved format state- 
ments must be simplified. 

Our future plans are as follows: 

1. Develop and implement an upward compatible 
trajectory oriented programming system. It must 
also be emphasized that all future programs written 
in the MSFC trajectory system will be FORTRAN IV 
compatible with other installations. 

2. Prepare abstracts for all modules both mathe- 
matical and program. 

3. Prepare a users manual for the system in- 
cluding complete description of statements and their 
source output. 

4. Provide training for the users. 

5. Study feasibility of adapting preprocessor to 
another machine. The Vectran Engineering Simulation 
System (VESS) preprocessor has already been trans- 
lated to ALGOL and will process FORTRAN IV state- 
ments on the Burroughs B-5500 computer. 

The trajectory oriented programming system 
planned by MSFC, Marshall Vectran Engineering 
Simulation System (MAR VESS), contains the only 
preprocessor which actually provides a system func- 
tion that will recognize a set of statements and create 
a trajectory program. No other system which we have 
studied can provide these necessary features. 


IV. AMTRAN 


I would like to mention at this time one other re- 
search effort which is underway at our Center and has 
as its goal the improvement of the man- machine re- 
lationship. This effort, which has been given the 
name Automatic Mathematical Translator (AMTRAN), 
is directed by Dr. Robert Seitz of the Research 
Projects Laboratory. AMTRAN is designed to be an 
automatic programming, on-line, multi -terminal 
computer system which should afford marked improve- 
ments in programming, debugging and turn-around 
times when it is fully developed. The system permits 
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a scientist or engineer to enter mathematical equa- 
tions in their natural mathematical format as they 
appear in a textbook and, barring complications, "to 
obtain an immediate graphical display of the solution 
on an output display device. The system is intended 
to be used for straightforward problem solution by 
the engineer or scientist with little computer experi- 
ence while at the same time providing the flexibility 
required by the experienced programmer to solve 
non-routine problems. A " sampler” version of the 
system is now available using a modified IBM 1620 
computer. 


! 

I V. RANDOM PROCESS THEORY 


The purpose of research in this area is to ex- 
amine the existing Computation Laboratory techniques 
used to reduce and analyze random process data 
toward the objective of devising new or improved 
applications of statistics and random process theory. 
The specific goals of this research are to reduce the 
data editing and computer usage time, to increase 
the "accuracy" of the statistical estimates of the 
processed data, and to recommend future applica- 
tions of existing data reduction equipment. These 
improvements are to be a result of the investigation 
of the techniques used by the Computation Laboratory 
and the appropriate application of: 

1. Digital filtering techniques 

2. Correlation function analysis 

3. Spectral smoothing techniques 
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The Computation Laboratory also initiated a re- 
search study with the Cornell Aeronautical Laboratory 
to do research in areas which would satisfy the 
Laboratory requirements of: 

1. Studying and applying the available random 
data processing techniques to the existing MSFC prob- 
lems, and 

2. Developing new and improved techniques of 
data processing. 

The following discussion indicates that the above 
requirements are being satisfied.* 

A. DIGITAL FILTERING 

Selection of an appropriate sampling interval 
which produces negligible frequency folding is para- 
mount to accurate digital data processing. The vast 
amount of literature available which describes digital 
simulation of transfer functions from the time re- 
sponse point of view can be used to produce pre- 
whitening filters having specific frequency character- 
istics. 

Taking the Tustin Transform of an analog notch 
filter will produce a digital filter which can be used 
for pre -whitening, with the possibility of total re- 
jection of one frequency. These notch filters contain 
relatively few weights. 

In situations where the power spectral density 
function of only a band of frequencies is of interest, 
digital heterodyning may provide a computational 
time savings in data processing. 


4. Special functions or processing 

5. Spectrum analysis of nonstationary functions. 

Research contracts were undertaken, to study 
numerical smoothing and differentiation methods. 

With these studies, digital filtering techniques were 
developed and investigated. The main effort was de- 
voted to linear digital (numerical) filters for per- 
forming smoothing, differentiation, and integration 
of discrete data and to do error analysis for these 
filters. 

The mathematical foundations were rigorously 
justified by beginning with classical Fourier theory 
and following through with the development of gen- 
eralized functions which led to specific functions used 
for filtering. This work is well documented in NASA 
Contractor Report CR-136. These desired digital 
filtering techniques were derived and are now being 
successfully applied to test data. 


B. CORRELATION FUNCTIONS 

After reading the analysis of different methods of 
estimating correlation functions, one should conclude 
that modifications should be made to any existing 
computational technique that does not consider both 
the accuracy of estimates and the computer time re- 
quired. Many types of correlation function estima- 
tors are given (autocorrelation being a special case 
of cross -correlation) . Extensive study of the "half- 
polarity" correlator is presented. Computer pro- 
grams are outlined, which will calculate, in minimum 
time, the "half-polarity" and "full -precision" cor- 
relation functions. It is also suggested that correla- 
tion computational techniques given in the reference 
are applicable. 


* Research Studies of Random Process Theory and 
Physical Application, NASA CR-61081. 
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C. OPTIMAL SMOOTHING OF POWER SPECTRAL 

DENSITIES (PSD) 

The appropriate application of proper techniques 
will produce spectral estimates with greater accuracy 
and also eliminate the need for pre-whitening of the 
signal prior to processing. 

In June, 1965, a 12-month extension to the 
project was initiated. The objectives are to extend 
and expand the techniques under the initial effort. 
Primary investigations will be the application of the 
non -stationary correlation function theory and digital 
correlation function computation techniques. The 
following list summarizes the technical effort and 
indicates the order of priority: 

1. Non- stationary data processing 

2. Stationary data processing 

3. Block diagrams covering application of data 
processing techniques developed. 

The following are the major accomplishments to 

date: 

1. The discrete data non- stationary correlation 
function theory has been developed. 

2, A solution for the form of the optimum filters 
to be used in the discrete data correlation function 
detector has been obtained. 


VI. DISCRETE OPTIMIZATION TECHNIQUES 


This laboratory has a contract with the University 
of Tennessee to do research in Discrete Optimization 
Techniques. The principal investigator is Dr. Gordon 
Sherman of the University's Computing Center and 
Mathematics Department. 

The problem is to maximize (minimize) a func- 
tion defined on a given finite set. Typical examples 
are: the shortest tour problem, the job shop sche- 
duling problem, and the transportation problem. 
Satisfactory solutions are available for some prob- 
lems of this class, while complete enumeration of all 
alternatives, if it were possible to do so, is the only 
known way of producing solutions for other cases. 

Dr, Sherman has taken a stochastic approach to 
the problem with the basic idea of combining intelli- 
gent search with random search. He has produced a 


family of algorithms that are quite efficient in the 
shortest tour type problem. These problems were 
used as test cases since the most research had al- 
ready been done on them. Detailed explanation of the 
method, algorithms, and results may be found in an 
article called "Discrete Optimizing" by Reiter and 
Sherman in the September 1965 issue of the Journal 
of Industrial and Applied Mathematics. 

VII. ANALOG COMPUTATION 
AND SIMULATION 

The traditional tool for simulation of dynamic 
systems has been the analog computer. The simi- 
larity between the real system and the program on 
the analog computer, and the possibility of identifying 
a block of the real system as a group of computer 
components, gives the simulation technique the ad- 
vantage of a model-like representation. This allows 
for an easy introduction of modifications and im- 
mediate observation of the effects of these changes. 
There are certain shortcomings in the use of analog 
computers. These are in the lack of random access 
memory, limited arithmetic precision, awkwardness 
in performing complex arithmetic, and others. These 
shortcomings led to a combination of the analog with 
the general purpose digital computers, thus preserv- 
ing the advantages of the analog while overcoming 
most of the shortcomings. This type of system is 
called a hybrid system. 

Hybrid computation, however, introduces prob- 
lems itself. Even though at many different places 
detailed investigations have been conducted, it was 
felt necessary to secure the support of an academic 
institution for basic studies in the area of error 
analysis of hybrid computation. 

Since this is a difficult and complex field, these 
studies are expected to become a long range effort. 
Some investigations have already been conducted by 
the Georgia Institute of Technology. The time limit 
for this review allows us to report only on the prob- 
lem area, the approach, and the more important 
results. Dr. Finn (of Georgia Tech) has investi- 
gated the errors introduced by sampling, by hold 
operation (zero and first order) for periodic , pulse 
shaped , and stationary band limited random functions 
of time. The sampling rate must be at least high 
enough to avoid fold over. When the highest frequency 
present in the continuous signal is f maximum, then 
the sampling frequency must be more than 2 f maxi- 
mum to avoid fold over. This is well known and fol- 
lows directly from a frequency presentation of the 
sampled signal. 
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When we concentrate on wide sense stationary 
random processes as time functions, we can make 
statistical predictions about the expected error and, 
what is more meaningful, the expected error square 
since the probability of positive and negative errors 
is equal. We can determine sampling rates, limiting 
the probability of our error to exceed a preset limit. 
The theory allows us to determine upper and lower 
bounds for the ratio mean square error for a band 
limited random process as a function of the sampling 
rate (zero hold) . These investigations are intended 
to be extended to higher order hold sampling tech- 
niques where similarly interesting results can be ex- 
pected. 

For sampling periodic functions, upper bounds 
for the ratio mean square error can also be given 
under the assumption that the initial phase is uni- 
formly distributed over all possible values in a ran- 
dom fashion. 

Dr. Finn has concentrated on investigating the 
error introduced by replacing the continuous function 
of time X(t) with a sampled representation using 
zero order and first order hold. Dr. Hammond, also 
of Georgia Tech, has worked on a system of n first 
order differential equations. With little loss in gen- 
erality, he starts with a class in which the first 
derivative is represented explicitly as a function of 
position and time. The hybrid computer uses its 
analog part for integration and its digital part for 
function generation. This allows one to derive for 
the error a system of linear differential equations. 
For short intervals the coefficients in these equations 
can be considered as constant and the forcing func- 
tion can be approximated by a staircase function. 

This not only allows one to use Laplace techniques 
for their analytical solution, but also provides a 
computer program of moderate complexity which can 
be incorporated with the hybrid program. A test pro- 
gram is presently set up in the Simulation Branch of 
the Computation Laboratory to study the usefulness of 
this approach. 


VIII. NUMERICAL INTEGRATION 


Because of the tremendous cost of modern com- 
puting equipment and the considerable amount of time 
used to perform certain types of studies, for example, 
orbit calculations, very substantial savings in com- 
puter time and dollars can be realized by even modest 
improvements in numerical integration techniques. 

The Computation Laboratory, in a continual search 


for better integration methods, has a dynamic pro- 
gram in mathematical and numerical analysis. This 
program is carried out by in-house staff members, 
specialists on a consulting basis, and through con- 
tracts with universities and some industrial firms. 


The numerical integration of differential equa- 
tions demands quite a large amount of computing 
time. Therefore, a great deal of attention has been 
given to devising more efficient methods of integra- 
tion. The laboratory has a research contract with 
.Vanderbilt University, Nashville, Tennessee, for 
the investigation of improved techniques for numerical 
integration of differential equations. The principal 
Investigator on the contract has been Professor E. B. 
thanks of the University's Mathematics Department. 
^Professor Shanks has devoted his efforts primarily 
'to a study of Runge-Kutta type processes. At the 
time the contract began, there existed the well-known 
Runge-Kutta formulas of fourth order requiring four 
evaluations of the function; the Kutta-Ny strom formu- 
las of fifth order requiring six evaluations; and the 
jless well-known Huta formulas of sixth order requir- 
ing eight evaluations. 

A paramount problem in trying to increase the 
order of the formulas in the Runge-Kutta sense is 
I that the number of conditions to be satisfied increase 
exponentially (essentially) and by the fact that the 
degree of the resulting algebraic conditions increases 
jby two at each stage; for example, a seventh-order 
jformula with nine evaluations involves 58 algebraic 
conditions with about half of them of twelfth degree. 

In such a complex system the notation becomes cum- 
bersome and a problem in itself. However, the prob- 
lem became tractable through use of the tensor calcu- 
lus notation. 

Dr. Shanks has been able to develop formulas of 
the sixth order with seven evaluations; seventh order 
|with nine evaluations; eighth order with twelve evalua- 
tions; and ninth order with seventeen evaluations. 

! 

By adopting a new view point in which not all of 
the algebraic conditions were exactly satisfied, Dr. 
Shanks has been able to develop formulas of fifth order 
accuracy with five evaluations; sixth order with six 
evaluations; seventh order with seven evaluations; 
and eighth order with ten evaluations. All experience 
to date indicates that these formulas are more ef- 
ficient than any of this type known previously. Addi- 
:tional details may be found in NASA Technical Note 
D-2920 and "Solution of Differential Equations by 
Evaluation of Functions," Mathematics of Computation , 
January 1966. 
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NEW ONE-STEP INTEGRATION METHODS OF HIGH ACCURACY 

By 

Erwin Fehlberg* 


SUMMARY ^ 

New numerical methods for the solution of 
periodic trajectories for the restricted three body 
problem are presented and factors affecting both the 
accuracy of results and the reduction of electronic 
computer time through the use of the new methods 
are discussed. Extension of the Runge-Kutta method 
to higher order of accuracy and the establishment of 
a pure series expansion method with transformation 
of the original differential equations to a second- 
degree algebraic system and application of recur- 
rence formulas has provided a method to more 
effectively use computer capability and point the way 
for use of the new methods in many space problems. 


I. INTRODUCTION 

The development of the electronic computer has 
created a need, and that need is becoming increasing- 
ly urgent, for more accurate, more powerful numeri- 
cal methods of computation. The computer is, of 
course, nothing more than a piece of hardware, how- 
ever complex. It obeys whatever numerical methods 
are programmed into it. Most numerical methods 
were developed long ago when the only tools available 
to mathematicians were pencil and paper, and perhaps 
a few tables of pre-calculated values (sines, cosines, 
logarithms, etc. ). The advent of the simple desk 
calculator helped. But even with the desk calculator, 
computational procedures had to be kept simple. 
Complicated operations, no matter how refined or 
necessary for the solution of certain complex prob- 
lems, were impractical or impossible. 

The rapid development of the modern electronic 
computer caught most mathematicians unprepared to 
use other than their old methods on the new hardware. 
Even now, the numerical methods used at many com- 
puter facilities are still the old desk calculator 
methods. For example, the standard method of 
Runge and Kutta, still widely used for the numerical 
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integration of differential equations, was developed 
around 1900. The method was quite suitable for use 
with hand-operated desk calculators. 

But many of the scientific and engineering prob- 
lems at Marshall Space Flight Center have become 
so involved that use of the standard, turn-of-the- 
century methods is completely out of the question. 

Not only are these methods often extremely slow, 
consuming excessive amounts of expensive computer 
running time; they are inaccurate, producing un- 
reliable results. Thus there is a pressing need at 
the Marshall Space Flight Center for new, more ad- 
vanced computational methods designed especially 
for use with the modern electronic computer. 

The Computation Laboratory has been actively 
seeking modern numerical methods suitable, in 
particular, for the solution of problems in astro- 
nautics and celestial mechanics. In recent years, 
several new approaches to the solution of ordinary 
differential equations have been developed in which 
such problems are expressed. One new approach 
is based on a power series expansion combined with 
a sophisticated, high-order Runge-Kutta procedure. 
Unlike the old methods still in widespread use, these 
new methods can conveniently be extended to any 
high-order accuracy desired. 

These powerful, high-order methods drastically 
reduce the errors involved in the numerical integra- 
tion of differential equations. Such errors originate 
both in the physical limitations of the computer, 
i. e. , round-off errors, and in the limitations of the 
numerical method programmed into the computer, 
i. e. , truncation errors. Moreover, in problems 
like the three-body problem, the new methods pro- 
ceed in large integration steps without impairing ac- 
curacy. Thus they are also much faster than conven- 
tional methods, which must proceed in extremely 
small steps to preserve some accuracy. 

The Marshall Space Flight Center is conducting 
extensive theoretical and numerical studies of 
periodic orbits of vehicles in the earth-moon system. 
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Data precisely defining a large number of such 
periodic orbits have been obtained using the new 
methods. Considering the effect of the truncation 
error, and to give some idea of the accuracy and 
speed possible with these methods, a sample periodic 
orbit in the restricted three-body problem was com- 
puted and found to retain its periodicity to within 0. 01 
millimeter (the distance from the earth to the moon 
is 384, 000 kilometers) . The computation took only 
about 5 percent as long as with the conventional 
Runge-Kutta-Nystrom method. 

These new methods can, of course, be used to 
solve many other problems in addition to problems 
in celestial mechanics. They are fully reported in 
the literature. 

II. AVAILABLE INTEGRATION METHODS 


A. MULTISTEP METHODS 

Methods for the numerical integration of dif- 
ferential equations are, broadly speaking, either 
multistep or one-step. Multistep methods were de- 
veloped as early as the nineteenth century, mainly 
for problems in astronomy. As their name indicates, 
these methods use the information from several back- 
ward computation steps ir calculating the solution for 
the current step. Muitistep methods (e.g. , the 
methods of Adams, Cowell, Gauss, etc. ) are very 
efficient for problems that can be integrated in steps 
of constant size. Since many such problems are 
encountered in astronomy, it is not at all surprising 
that a number of multistep methods have been de- 
veloped by astronomers. Multistep methods also have 
the very great advantage that they generally require 
only one or two evaluations of the differential equa- 
tions per step, and they can be extended to any order 
of accuracy simply by adding higher-order difference 
terms to the formulas. Hence, unlike the Runge- 
Kutta method, which is the classical one-step method, 
they are quite fast on an electronic computer; they 
are economical, and they can be made very accurate. 

But multistep methods do have a number of 
major disadvantages. They are not self-starting, 
but require a special starting procedure. A history 
of known values is needed before computation can 
begin. Thus a number of backward values must be 
created by means, for example, of an iterative pro- 
cedure. And if the integration step size has to be 
changed during the computation, if, for instance, the 
step size must be reduced to preserve accuracy, 
additional time-consuming iterations are needed to 
build a new difference scheme— time-consuming 


because the iterated values for the changed step size 
must be of high accuracy. For these reasons, among 
others, multistep methods are largely restricted to 
problems that can be integrated entirely in steps of 
the same size. This is, of course, the case in the 
determination of astronomical orbits, where the 
distances between the attracting bodies are changing, 
but not radically. 

This is not at all the case in, say, the interest- 
ing orbits of the restricted three-body problem. 
Figure 1 shows a typical periodic orbit of the re- 
stricted three-body problem. The earth and the 
moon are the two attracting masses. These and the 
space vehicle are shown in the rotating coordinate 
system, in which the x-axis always extends from 
the earth to the moon. Every fifth integration step 
is indicated, except in the vicinity of the earth, 
where the steps were too numerous to show. Ob- 
viously, the integration steps in the vicinity of both 
the earth and the moon are much smaller than those 
where the vehicle is far from either attracting body. 


t*4 



FIGURE 1. PERIODIC ORBIT OF THE 
RESTRICTED THREE-BODY PROBLEM 

This kind of flexibility— to be able to increase the 
step size as much as possible or to be able to de- 
crease it as much as necessary— is essential for 
efficient integration of such problems as the re- 
stricted three-body problem. It speeds the integra- 
tion during the large part of the orbit where the space 
vehicle is near one or the other of the attracting 
masses. 

This need for flexibility in the size of the inte- 
gration step is caused by two kinds of error that 
accumulate during a computation. Too large a step 
size results in an unacceptably large truncation er- 
ror because truncation error is proportional to a 
certain power of the step size. Too small a step 
size not only slows the computation, thus wasting 
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second degree. For special differential equations 
the procedure has been outlined in earlier papers by 
Steffensen [ 4] , Rabe [ 5] , and the author [ 3] . The 
procedure is based on the fact that the high-order 
derivatives of a second-degree system of differential 
equations can be conveniently obtained on a computer 
by recurrence formulas. 

This can best be illustrated by the transforma- 
tion procedure in a simple example. Consider the 
differential equation 


Introduce the auxiliary function 
-x 

e = u (2) 


at At can be adjusted immediately in such a way that 
| X n +i (At) 1 * 1 1 remains within prescribed limits. 
(For safety it might sometimes be advisable to con- 
sider more than just one term of the truncation 
error. ) Unlike Runge-Kutta or multistep methods, 
no repetition of any computation is necessary if the 
step size fails to meet the requirements for the 
magnitude of the truncation error. The Computation 
Laboratory knows of no other method that offers such 
easy step-size control. 

In this simple example, there is no real need to 
introduce auxiliary functions, since a repeated dif- 
ferentiation of the differential equation ( 1 ) can be 
performed without difficulty. A more representative 
example follows to illustrate how convenient the 
method can be. 


and obtain from equations (1) and (2) a system of 
second-degree algebraic differential eaquations 

(3) 

Substituting the power series expansions 


dx 

dt 


= u. 


du 

dt 


= -u 2 


IV. POWER SERIES EXPANSION METHOD 
APPLIED TO THE RESTRICTED THREE-BODY 
PROBLEM 

Clearly, the following equations, for the re- 
stricted three-body problem in the rotating coordinate 
system, are not nearly so simple as equation (1): 


x= E x „ • (t-to)". u = £ U • (t-to) V (4) 

y=0 


into equation (3) and comparing coefficients for the 
terms with (t - t 0 ) n results in the following recur- 
rence formulas for the coefficients in equation (4) : 


(n+1) X . = U 

' rtf-1 n 


(1*1) U = - T U • U 

n+1 Li v n - v 
v=0 


(n=l ,2,3,...) 

(5) 


Since the first coefficient Xq is known from the 
initial value x(t c ) for the step and the first coefficient 
U Q can be obtained from equation (2) , all following 
coefficients X^, (v = 1, 2, 3, . . . ) can easily be 
computed from the recurrence formulas of equation 
(5) , a very convenient procedure for electronic 
computers. 


d 2 x , dy , , v x + u 

dt 2 " X 2 dt ( ' [(xt-jutf + y 2 ] 3/2 

_ X - (1 - u) 

(x-i+|U) 2 + y 2 ] /2 

dV „ dx , v 

dt 2 y 2 dt ^ [(x+jn) 2 + y 2 ] 3 / 2 

" M [( x-l+J) 2 +y 2 ] 3 / 2 

where p = the relative mass of the moon in the earth- 
moon system. 

There exists a first integral of these equations 
of motion, the so-called Jacobi integral 

J = i [(sf) 2+ (dt - ) 2 " x2 ~ y2 ] ~ [f x-1+^7 2 +y2"] ^ 2 " 



It is quite obvious, too, that the power series 
expansion method allows for an extremely simple 
automatic step-size control. Assuming truncation 
of the expansion in equation (4) for x after the term 
Xn(t-t 0 ) n , leading term of the truncation error 
of x can easily be found by extending the computation 
to the next coefficient X^j. If the truncation error 
turns out to be too large or too small, the step size 


— ^ v 9 i / 9 — * Const. 

I (x-l+g ) 2 + y 2 ] V* 

Auxiliary functions are again introduced 
r 2 = (x+m) 2 + y 2 , s 2 = (x-l+ft) 2 + y 2 ) 


(7) 


( 8 ) 
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expensive machine time, even more seriously, it 
results in an unacceptably large round-off error be- 
cause round-off error is a direct function of the 
number of steps that must be taken to compute the 
entire problem. 

For use with computers, the ability to selectively 
vary or to automatically control the integration step 
size is fundamental to the efficient, fa st-and- accurate 
integration of problems like the restricted three-body 
problem. This cannot be done with multistep methods 
except at considerable expense in increased complex- 
ity and increased computer running time. 

B. ONE-STEP METHODS 

One -step methods lend themselves more readily 
to step-size variation. In fact, the step size can be 
changed at any time and can immediately be accom- 
modated to local conditions at any point in an inte- 
gration, One-step methods are also self-starting. 

But they, too, have a number of disadvantages. The 
classical Runge-Kutta method is of only fourth-order 
accuracy. Several extensions have been made in the 
last decade, notably by Shanks [ i , 2} , but Shanks 1 
method, the most accurate Runge-Kutta procedure 
developed to date, is still of relatively low-order 
accuracy (up to eighth-order) . Runge-Kutta methods 
also tend to be slow because of their great complexity. 
The differential equation must be evaluated many 
times for each integration step. Shanks 1 eighth-order 
method, for example, requires 12 evaluations per 
step. If, in addition, the differential equations are 
complicated, if they contain transcendental functions 
(sine, cosine, exponential functions, etc.), the 
method becomes excessively slow and the cost in 
computer running time will certainly be high. 

Another weakness of Runge-Kutta methods is that 
they, too, like the multi-step methods, lack an 
economical procedure for automatically adjusting the 
step size to the local conditions of the problem. As 
with all one-step methods, the step size can be 
changed at any time, but no economical control pro- 
cedure seems to exist to do it automatically. In fact, 
using Runge-Kutta methods, one never knows whether 
the proper step size (for a combination of maximum 
accuracy and speed) is being used. There is no easy 
way to determine this. Apart from somewhat doubt- 
ful rule-of-thumb control procedures, there exists 
only Richardson’s well-known method of the deferred 
approach to the limit. A step is computed, recom- 
puted with half the step size (or double the step size), 
and then, by an extrapolation procedure, the results 
of the two computations are compared. This, how- 
ever, doubles the computational effort or, more 
exactly, doubles computer running time merely for 
the benefit of step-size control. 


This, briefly, was the state-of-the-art in 
methods for integrating differential equations when 
the Computation Laboratory began research in the 
field. The available methods were rather inaccurate 
and slow, and costly in machine time even for the 
solution of a relatively simple problem like the re- 
stricted three-body problem. For more complex 
problems in orbital mechanics, like the n-body prob- 
lem, they could well turn out to be prohibitively slow 
and, worse, intolerably inaccurate. An original, start- 
ing 16-digit accuracy could easily dwindle to two or 
one or no accurate digits at the end of a long, com- 
plex computation. Also, it seemed absurd to pay 
heavily for auxiliary features like step-size control. 

III. POWER SERIES EXPANSION METHOD 


As a first quick improvement, the range of prob- 
lems normally solved by pure power series expansions 
was expanded [ 3] . The use of pure power series ex- 
pansions to solve differential equations is not exactly 
new, of course.* But unless the differential equation 
under consideration was extremely simple, pure 
power series expansion methods had generally been 
discarded as leading to cumbersome and lengthy com- 
putations for the derivatives they require. And, in 
fact, the method suggested requires a repeated total 
differentiation of the differential equation (s) with re- 
spect to the independent variable to obtain the neces- 
sary coefficients of the power series expansion. 

Only a few years ago, the repeated total different- 
iation of a differential equation was not considered 
feasible, since, with increasing order, the deriva- 
tives become rather unwieldy expressions. But with 
fast electronic computers such a procedure is gener- 
ally quite feasible. It is well-known that in the last 
few years considerable progress has been made in 
the automatic differentiation of formulas by computers. 

Moreover, as an even more effective approach, 
apart from a straight-forward differentiation of the 
differential equation (s) , a great number of these can 
be differentiated in a rather simple way by first trans- 
forming them, through introduction of auxiliary func- 
tions, into algebraic differential equations of the 


* All commonly used integration procedures are, in 
fact, based in part on a power series expansion. The 
coefficients of Runge-Kutta formulas or the coeffic- 
ients of multistep methods are obtained by expanding 
in power series (Taylor series, etc. ) and then find- 
ing systems of equations of condition for these co- 
efficients. 


11 





Introducing equation (8) into equation (6) trans- 
forms the original system into the following second- 
degree algebraic system, which can be integrated 
directly by power series expansions: 


d^x dy 

W = x + 2 dt " u(x+/i) " v ( x_1+ ^) - 

d 2 y dx 

dt? = y - 2 ^ - vy 


du 0 dr _ dv , _ ds 

r— +3u — = 0, s — + 3v — = 0, 
dt dt dt dt 


r 2 - (x+ij .) 2 + y 2 , s 2 = (x-l+pi) 2 + y 2 



Again, of course, new differential equations for 
the auxiliary functions have been added, but the new 
system is completely algebraic, containing only pro- 
ducts of two functions throughout. 


Thus all the Taylor coefficients of X and Y can be 
obtained. These expressions can be extended to as 
many terms as desired. There is no restriction 
whatsoever on the order of the formulas. This is in 
distinct contrast to Runge-Kutta formulas, in which 
each advance of only one order in accuracy has taken 
many years to establish and, as mentioned earlier, 
the highest known order is only the eighth. There is 
no such problem with recurrence formulas. It is such 
recurrence formulas, evaluated automatically on the 
computer and extended to any order desired, that form 
the basis of both the pure power series expansion 
method and, as will be shown in the next section, the 
improved Runge-Kutta method developed by the 
Computation Laboratory. By a suitable transforma- 
tion the original differential equation(s) is reduced 
to a second-degree algebraic system and then the re- 
currence formulas are applied whose coefficients are 
determined automatically on the computer. 


Let the power series expansions be 
x=V X (t-t)", y=T Y (f-t )" \ 

v =0 v ° i /=0 V ° J 

U= Z U (t-t 0 ^ v=£ V (t -t/ \ (10) 

y=0 y=0 I 

r = y R (t-t )", s = y S (t-t ) V / 

u n V o' v 'o' / 

v ~ 0 v =0 1 


The Computation Laboratory has solved many 
problems— restricted three-body, motion of an elec- 
tron in the field of a magnetic dipole, and others— by 
a pure power series expansion method. However, 
while the method proved superior to other existing 
methods, there still seemed to be some room for 
improvement. For example, the number of terms in 
all the sums in equation (il) increases with increas- 
ing n. Hence computer running time gets longer for 
higher-order coefficients. 


The first coefficients X 0 , X t and Y Q , Yj, are 
known at the beginning of the integration step. The 
forst coefficients R Q , S 0 , U 0 , and V Q are then de- 
termined from equation (8). 

Inserting power series expansions from equation 
(10) into the original system of equation (9) , the 
following recurrence formulas for the succeeding 
coefficients are obtained: 


“ 11 n-i 

28 B - ) X X + 2(xX +V Y Y - Y R R 

° n v n ~ u n "o * n_t; „=1 

n n n-1 

2S S = A X X - 2 ( l-/jt)X + y Y Y -V ss 

° n ^0 n Jb » n-v ™ 

n n-1 

nR u „ = -3 V i/R U -V i-u R 

° n v n-y v n-y 

n n-1 

nS v = -3 y 1 /S V - V yV S 

° n u n -" J = i V 

(n+1) nX nM = X n-l + 2 nY „-‘ ,U „-l + < 1 -'') V „-l Y 

i’~Q 

n-1 

(n+ 1) nY = Y - 2nX - Y <U + V )Y . 

n+ 1 n-1 n •— < y y n-l-y 


\ 


(n = 1,2,3,...) 

\ (H) 


It is well-known, too, that power series expan- 
sions have certain limitations with respect to the 
truncation error. When one truncates the expansion, 
there is no way of covering the remainder of the error, 
which is roughly equal to the leading term. This is 
unavoidable in a power series expansion, although it 
is even more of a problem in multistep methods. This 
need not be such a problem with Runge-Kutta methods, 
and it is for this reason that the Computation Labora- 
tory has developed the Runge-Kutta transformation 
method, which combines the high-order accuracy of 
power series expansions with a good coverage of the 
truncation error. In fact, since, to a certain extent 
at least, the leading term of the truncation error can 
be covered, the Runge-Kutta transformation method 
radically reduces the truncation error. 


Thus a combination of the two methods should 
not only be more accurate than the pure power series 
method, it should also provide an advantage in speed 
because larger step sizes can be used. The combined 
method is described in the next section. 
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V. RUNGE-KUTTA TRANSFORMATION 
METHOD 

A. FORMULAS OF ANY DESIRED HIGH ORDER 


Consider only second-order differential equations 
since it is these that are most frequently encountered 
in physics and mechanics. (The method works for 
first-order systems as well [ 6] . ) For brevity, 
formulas for only a single equation will be written, 
although the method holds in exactly the same way 
for systems. Letting x be the original dependent 
variable, then: 


x = f(t, x, x) 


x(t 0 ) = x 0 , x(t 0 ) = x 0 


( 12 ) 


Next a transformed variable x T is introduced, 
which equals the original variable minus the first 
m+2 terms of the power series expansion for x 

m+2 

X T = x " Z X „ (t - t.)' ( 13 > 

v=i 

m+2 

x T = X - Yj v \ (t - to)- • (14) 

V=1 


f T (t 0 + 

ajh, 

x 0 , 

0) h 

f T (t 0 + 

a 2 h, 

*0 

+ /3okjh 

f T (to + 

«3 h > 

*0 

+ yokjh 

6jk 2 )h 




f T (t 0 + 

<*4 h > 

x o 

+ eokjh 

J 1 k 2 + i 7l k 3 )h 




(16) 


and 

x T = x 0 + (Cjkj + C 2 k 2 + C 3 k 3 ) h + 0 (h m+ J ) 

x T = 0 + Cjkj + C 2 k 2 + C 3 k 3 + Ofh™* 5 ) 

x T = x 0 + (Cjkj + C 2 k 2 + 6 3 k 3 + 6 4 k 4 )h + Ofh”^ 6 ) 


(17) 

Three substitutions yield an accuracy in x, x to the 
^m+4th ^ erm> w here m is the number of differentia- 
tions performed in equation (13) before x can be re- 
placed by xj. The fourth substitution in equation 
(16) yields the truncation error term required for 
automatic step-size control. 


Performing this subtraction results in a function 
with zero derivatives for t=t Q up to the m+2nd order. 
The following differential equation is obtained for the 
transformed function x^: 

m+2 \ 

k T =f T =f “X V ( l ’ _1 ) Xj, (t - to)^ / 


x T (t 0 ) = x (t 0 ) = x 0 , x T (t 0 ) = 0 


That the first m+2 derivatives equal zero con- 
siderably faciliates establishment of Runge-Kutta 
equations of condition. Furthermore, accuracy 
can be of any high order desired simply by subtract- 
ing enough terms from the original function. In 
other words, a very simple function is always created 
that has zero derivatives up to the m+2nd order. 
Runge-Kutta formulas of any high order desired can 
then be obtained for this transformed function merely 
by choosing m sufficiently large. For example, the 
Runge-Kutta formulas for the transformed differential 
equation (15) just given would read 


Thus an additional advantage of this approach, 
which also distinguishes it from any other Runge- 
Kutta formulas, is the very simple, economical (in 
computer running time) procedure for control of the 
truncation error. Only the first three evaluations 
are needed for the actual computation; the fourth 
evaluation gives an improved value for k which is 
accurate to one further power of h. By subtracting 
these terms, the leading term of the truncation error 
is represented with sufficient accuracy 

T x n T -x T =[ (C, - C^kj + (C 2 - 6 2 )k 2 

+ (C 3 - 6 3 )k 3 - 6 4 k 4 ]h . (18) 

Full details on these new high-order Runge- 
Kutta formulas are given in reference 6. 

B. FORMULAS WITH AN ARBITRARILY SMALL 

TRUNCATION ERROR 

It may be noted, without going into detail, that 
in more recent work the Computation Laboratory has 
established Rungf -Kutta formulas in which a param- 
eter a and the absolute value of all members of the 


14 



ertin fehlrerg 


leading term of the truncation error, for x^ as well 
as for Xrp, can be made as small as desired, but not 
zero since some coefficients would then become 
infinite. Full details on these high-order formulas 
with an arbitrarily small truncation error are given 
in reference 7. 

These new formulas required again a suitable 
transformation of the original differential equations. 
This transformation is based on a power series 
expansion. Any desired degree of accuracy can be 
obtained by doing only three or four evaluations of 
the differential equations. This is not possible with 
other Runge-Kutta type formulas now found in the 
literature. There is a small penalty to be paid in 
computer time since the recurrence formulas must be 
evaluated. The additional computation is not great 
and, because the method is of high order, the inte- 
gration can proceed in larger steps without impairing 
accuracy. The computer running time is much faster 
than for other known integration methods. This is 
shown in Section VII. 


VI. SOME OTHER MODERN RUNGE-KUTTA 
FORMULAS 


Briefly, for comparison, consider two other 
modem Runge-Kutta methods: the Shanks explicit 
method and the Butcher implicit method. To simplify 
comparison, both methods are presented in eighth- 
order form. As mentioned earlier, Shanks’ formulas 
are available only to the eighth order. Butcher’s 
implicit formulas are available to any order. 


First, consider Shanks ! explicit formulas [1,2] 
for x = f (x) \ 

ki=f(x 0 )h 1 

k 2 = f(x 0 + « 2 iki)h I 


k 3 = f(X 0 + Q!31 k l + a 32 k 2) h 


k 4 = f (x 0 + a 4 i k i + a i2 k 2 + a 43 k 3 )h 


k i2 f( x o + a 12> l k l + • • • + a i2> ll k ll) k 


12 

x = x o + Z C „V 0 ( h 9 ) 

1 


(19) 


Each integration step here requires 12 substitutions, 
kj through k 12 , which are multiplied by certain weight 
factors and summed to obtain the new value for x. 

But because these formulas include no procedure for 
controlling the truncation error, each integration 
step really requires 23 substitutions if Richardson’s 
extrapolation procedure is used for step-size control, 
i. e. , 2- 12 substitutions with one substitution omitted 
since the first substitution occurs twice in the com- 
putation. 


In 1964, Butcher [8, 9] published two noteworthy 
papers on implicit Runge-Kutta methods. Following 
are his eighth-order formulas for x = f(x): 


k i = f ( x 0 ) h . 

,k 2 = f(x fl + /?21 k l + 0 22 k 2 + /?23 k 3 + 024^4) k 

k 3 = f( x 0 + ftu k l + As2 k 2 + 033 k 3 + /3s4 k 4 ) h 

k 4 = f( x 0 + 0 41 k i + 042 k 2 + 04S k 3 + ft 4 k 4 )h \ . (20) 

k 5 = f ( x o + &l k l + 052 k 2 + 053 k 3 + /3 M k 4 ) h \ 

x = x o + Z c i k , + °( h9 ) / 


Unlike explicit formulas, where the increments 
for x: kj, k 2 , k 3 , etc. , are successively computed, 
sach value depending only a previous values, implicit 
formulas require an iterative computation. Any 
increment kp depends not only on the preceding incre- 
ments k lf k 2 , ... k v _i but also on k^ itself and on the 
succeeding increments k^+i, k v + 2 , .... Naturally, 
this iterative computation is more involved than the 
straightforward procedure for explicit Runge-Kutta 
formulas. But implicit formulas do require con- 
siderably fewer substitutions than explicit formulas. 
Formulas (20) require only five substitutions, only 
three of which are iterative, per step versus 12 for 
the comparable Shanks’ formulas. However, the 
iteration tends to be slow. This is demonstrated in 
Table I. 
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TABLE I 

COMPARISON OF TWELFTH-ORDER METHODS, 
RESTRICTED THREE -BODY PROBLEM. ( 


Method 

Final x 

Final y 

Number 
of Steps 

Computer Running 
Time (min) 






rkb (2) 

1.20000 00000 00013 

-1. 04935 75098 30328 

216 

0.88 

pse (3) 

1. 19999 99999 99981 

-1. 04935 75098 30303 

493 

0. 21 

RKT 1 (4) 

1.20000 00000 00001 

-1. 04935 75098 30321 

389 

0. 15 

RKT 2 (5) 

1. 20000 00000 00013 

-1. 04935 75098 30332 

290 

0. 13 


(1) Corresponding results for Shanks 1 eighth-order formulas are: 1.20000 00000 00002; - 1.04935 75098 30310; 
814 steps; and 0.46 minute. Note that eighth-order formulas are, of course, not competitive in speed or 
accuracy with twelfth-order formulas. This is far more obvious in more complex differential equations, 

^s can already be seen in the more complex, but still relatively simple, case of the restricted four-body 


problem. 

(2) RKB = Runge-Kutta-Butcher method [ 8, 9] 

(3) PSE = Power series expansion method [3] 

(4) RKT 1 = Runge-Kutta-transformation method [ 6] 

( 5) RKT 2 = Runge-Kutta-transformation method [ 7] 


VII. CONCLUSIONS 


As an example of the computing speeds possible 
with the modern methods described, the periodic 
orbit shown in Figure 1 was computed. The orbit 
has the following initial values: 

x 0 = 1. 2, y 0 =0, x 0 = 0, y 0 = -1. 04935 75098 30320 

(M= 1/82.45) . 

( 21 ) 


To preserve sufficient accuracy for a true compari- 
son of the methods, the initial value y Q was computed 
in 20-digit arithmetic. The computations were exe- 
cuted on an IBM 7094, Model II computer (16 digits) . 
Table I shows the results for one complete orbit. 

For all methods compared in Table I, we lose 
about two digits on a 16-digit computer. The Butcher 
method, which uses the fewest integration steps, is 
extremely accurate but it is also extremely slow. 

In all cases, however, the deviations are negligible, 
being of the order of 0. 01 millimeter for this particu- 
lar orbit. But the method is nearly seven times as 


fast as the Butcher method, i. e. , it would cost about 
seven time as much in computer rentals to use the 
Butcher method. Thus, even in a relatively simple 
problem like this example, it pays to use the most 
efficient integration method available. 

A group, headed by Mr. Mert C. Davidson, is 
being set up in the Computation Laboratory to explore, 
in detail, applications of these methods to practical 
problems that will exploit their full possibilities. For 
example, a program has already been written to solve 
the complete n-body problem (including oblateness 
terms) , as a whole, with no reliance on data from 
relatively inaccurate external sources like ephemeris 
tables, etc. The computation of this problem has been 
rapid and highly successful. But many new worth- 
while applications still need to be developed and ex- 
ploited. 

Finally, it should be noted that the subject matter 
of this paper is given much more thorough coverage 
in the paper New One-Step Integration Methods of 
High-Order Accuracy Applied to Some Problems in 
Celestial Mechanics , which will be published shortly 
by NASA. 
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CELESTIAL MECHANICS 


SUMMARY 


2 93*7 


For the restricted problem of three bodies the 
perivation of periodic solutions different from the 
classically known ones is discussed and the ideas 
used for their existence proofs are briefly outlined 
and related to classical methods. In particular, the 
derivation of closed perturbed precessing elliptic 
orbits of arbitrary eccentricity and small major axis 
about the smaller one of the two attracting bodies 
with arbitrary mass ratio is indicated. Two numer- 
ical examples of recently discovered closed trajec- 
tories are included. 


I. INTRODUCTION 


Among the problems of celestial mechanics, 
which are important for space flight applications, the 
restricted three body problem plays a central role. 
This problem is concerned with the description of 
the possible trajectories of a particle (manned or un- 
manned satellite, meteorite, planetoid) of negligible 
mass under the gravitational attraction from two 
heavy celestial bodies, which are assumed revolving 
according to Kepler's laws on circles about each 
other. Limiting our attention to the two-dimensional 
case the equations of motion of the particle can be 
written in the form 


x + 2ix - x = -^( x + ^i) |x + ^r 3 -pi(x- v)\x-v\ 3 , 


Hill, Poincarg, Birkhoff and others, and is still 
being pursued vigorously. It will be our sole con- 
cern in this presentation. The study of periodic 
solutions of equation (1) is of interest for several 
reasons. First, since the restricted three body 
problem presents a non-integrable dynamical system, 
every contribution toward an understanding and a 
description of its general solution afforded by par- 
ticular solutions is highly welcome. Second, recent 
advances by Kolmogorov, Moser and Arnold in the 
areas of stability and almost periodic motions have 
led to an understanding of the behavior of dynamical 
systems in the vicinity of its periodic motions. And 
third, some periodic solutions of equation (1) are of 
great practical interest in dynamical astronomy or 
in space flight mechanics. 

Our present knowledge of periodic solutions of 
equation ( 1) is still modest. Without discussing the 
classically known solutions, which are either near 
the libration points, or are close to circular solutions 
(for small \x > 0, or, for arbitrary ji, when near one 
of the masses or far away from both masses) , or 
which are inside a closed zero-velocity oval about 
the heavier mass closing only after many revolutions, 
etc. , we will give a description of some recently 
discovered periodic solutions and of the ideas used 
for their existence proofs. These new solutions are 
characterized by their relationship to Keplerian 
elliptic motions of positive and possibly large eccen- 
tricities, presenting relative to equation (1) , a 
situation which classical researchers attempted in 
vain to illuminate although they had essentially 
created the methods with which to attack such pro- 
blems. 


C = d/dt , i 2 = -l) (1) 

where x = x A + ix 2 is the complex position vector of 
the particle, referred to a rectangular coordinate 
system, which rotates with unit angular velocity 
about the center of mass of the two heavy bodies with 
masses n and v = 1 as origin. 


II. PROBLEM 

Let us describe our problem. Equation ( 1) 
approximates the equation of motion for the Kepler 
problem (two body problem with one mass trans- 
formed to rest) 


The attempt, among others, to exhibit periodic 

solutions of equation (1) has received great empha- z = -mz|z|~ 3 , m > 0 (2) 

sis and led to some success through the work of 


* Staff Scientist, Computation Lah ; . MSFC. 


in an inertial coordinate system after a rotation 
z = e^x ; for example , if (x > 0 is very small, m = v 


19 


! 


RICHARD F. ARENSTORF 


and the particle does not approach the smaller body 
too close (planetary case) , or when the particle 
moves in the vicinity of the smaller one of the attrac- 
ting bodies ( satellite case) , etc. Now equation ( 2) 
has elliptic solutions, say with major half axis a > 0 
and eccentricity e, 0 < e < 1. Such a solution is 
closed in the rotating x 1# x 2 coordinate system, if 
its period T 0 = 27ra 3 / 2 is a rational multiple of the 
period of e^, i. e. , if a = (m/k) 2 / 3 with integers 
m, k ^ 0. Under this assumption the resulting ro- 
tating elliptic orbit to be described by x = x v ( t) is 
closed after k - m revolutions about the origin, having 
the period T* = 27rm. The problem is to find periodic 
solutions x = x(t) of equation ( 1) which are near 
x* (t). 

This problem can, in the simpler planetary case, 
be solved with the classical methods devised by 
Poincare. The only additional idea needed consists 
of the application of an appropriate periodicity con- 
dition instead of the classical condition of return of 
the motion to its initial state after time T > 0. This 
classical condition leads to a singular case, even 
after reduction with the help of the Jacobian integral 
of equation ( 1) , when applied to the generation of 
periodic solutions of equation (1) from the above 
x* (t). The difficulty can be overcome by using the 
condition 


x(t) = real, x(t) = pure imaginary at t = 0 
and t = ±T > 0 


(3) 


for a solution x = x(t) of equation ( 1) . Equation ( 3) 
implies that the curve x = x( t) (0 < t < T) becomes 
symmetric about the real axis of the x-plane and 
closed. Thus, x(t) becomes periodic with period T. 
To satisfy equation ( 3) , the solutions x of equation 
( 1) represent analytic functions not only of t, but 
also of the parameter n in equation ( 1) and of the 
initial position and velocity coordinates 


x.(0) = £ , x . ( 0) = t). , (j = 1,2) . (4) 

J J J J 

Then equation ( 3) can be rewritten as 


x 2^ £ 1» $2» Vl> T?2»M) - x i( ^T, £ 2 , 772,M) - 0, 

£2 = m = 0 (5) 


giving two scalar real equations for the unknowns 
T, £*, rfe. When fj. = 0, these equations have a known 
solution (say T* , , 172*) belonging to the genera- 

ting solution x* (t) of equation ( 1) with ju = 0 after 
proper choice of its initial values. Since the re- 
spective Jacobian determinant with respect to T and 
772 does not vanish (to establish this fact constitutes 
the decisive part of the existence proof) , the im- 
plicit function theorem leads to the existence of 
solutions for T, 77 2 of equation (5) near T* , , 

772* for sufficiently small ijl > 0 and thus to periodic 
solutions of equation (1) near x* (t) [1]. 

By proper choice of m , k, and e above the pro- 
cessing elliptic orbit (t), (0 < t ^ T' fi ) can be made 
to pass the attractive bodies at prescribed small 
distances and this property will still hold for the 
resulting periodic trajectories x(t), (0 < t < T) of 
equation (1) with fx > 0, since y is small. Such 
trajectories are of great astronautical interest for 
space flight in the Earth-Moon system. 


III. SATELLITE CASE 


We now come to the more difficult satellite case 
of our problem. This time we need new ideas to show 
that periodic solutions x(t) of equation (1) near 
x*(t) exist (since 1 > 0 is small, for instance) and 
the periodic motion is to take place in the near 
vicinity of the body of mass recalled planet) , where 
the disturbance exerted by the other more massive 
body, being nearly at rest in the inertial coordinate 
system , causes large deviations from Keplerian 
motion for the third body near the revolving planet. 

If x* (t) is simplified to a circular solution by putting 
c = 0 and dispensing with the condition that Tq/2tt be 
rational, our problem has been solved already in 
different ways by Hill, Brown, Moulton, Wintner and 
Siegel. But these results, which are based on power 
series expansions of the coordinates xj, x 2 in powers 
of the small period, give no indication of the exis- 
tence of periodic solutions near x*(t) with e > 0 and 
small a = (m/k) 2 / 3 . The first result in the direction 
of the present problem, although only in the planetary 
case, was obtained by Birkhoff using the Poincare - 
Birkhoff fixed point theorem for annulus mappings, 
and more recently by Moser, who used the Birkhoff - 
Siegel fixed point theorem for local area-preserving 
mappings to get a more accurate description of the 
location of the obtained solutions of equation (1). A 
similar result in the satellite case, giving for each 
sufficiently small value of the Jacobian integral the 
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existence of countably many periodic solutions which 
close only after many revolutions about the planet, 
was recently established by Conley using the Poincare- 
Birkhoff fixed point theorem and a new construction of 
the classically known nearly circular solutions men- 
tioned above. These periodic solutions of equation ( 1) 
still lack a more accurate geometrical description. 
Adequate references to the literature are contained 
in references 1 and 2. 

We shall now give a brief description of the ideas 
and techniques which lead to the existence of countably 
many families of 1 -parametric solutions x(t) near 
x* (t) in the satellite case of equation (1) , belonging 
to sufficiently small values of m/k - a 3 / 2 and having 
the family parameter e which ranges over suitable 
closed intervals contained on 0 < e < 1. These solu- 
tions exist for all 0 < v < 1. 

We transform equation (i) by a translation w = 
x + fji into 

w + 2iw - w + ^ w | w | 3 

= -p(l + (w - 1) | w - 1| -3 ) =P q (w). (6) 

Here the right hand term P Q (w) is the disturbing 
function, which vanishes at w = 0, i. e. , at the location 
of the small planet of mass v m At w = 0, the left side 
is singular, however. Replacing Po(w) by 0 in equa- 
tion (6) leads, after the rotation z = eitw, to equation 
(2) with m = v. Thus equation (6) is close to the 
integrable Kepler problem for small |w| even though 
[x = 1 - v is not small. Again, use of the periodicity 
conditions of equation (3) , with x replaced by w, is 
decisive and assures a non-vanishing Jacobian rela- 
tive to the unperturbed elliptic motion x* (t). Re- 
writing equation ( 3) in the form of equation ( 5) , 
however, is of no use since now fx is not considered 
as a small available parameter, but is fixed and 
nearly 1. Despite this we shall solve equation (3) 
with equation ( 4) for T and 772 * as in the case of 
equation (5) , with the help of an implicit function 
theorem by application of the following idea. 

We replace in the right side of equation (6) the 
given function P Q (w) by an arbitrary function P(w) 
from a suitable set F of functions, which contains 
P 0 (w) especially. The resulting solutions w( t) (or 
x(t) = w(t) -ix just as well) then depend upon their 
initial values and upon P(w). Therefore equation (3) 
can be rewritten in the form of equation (5) , but with 
M replaced by P if P is the name of the chosen function 
P(w) from F. Now P can be considered as a gener- 
alized parameter varying over F, instead of the real 


parameter 11 . When P = 0 (the zero element in F) , 
equation (5) again has a known solution (say T* , , 

tj* 2 ) determined by x* (t). Thus, using an appropri- 
ate implicit function theorem we arrive at the exis- 
tence of solutions T, , n 2 of equation (5) near T* , 

rf ' 2 f° r sufficiently small P of F (say for || P j| < 
r*) after having introduced a suitable norm ]| . . || on 
F. But then it is decisive that the given P Q in equation 
(6) satisfies ||P 0 || ^ r * to obtain periodic solutions 
of equation (6). Since P Q is not available to choice, 
a lower estimate is needed for r* and not merely the 
existence of an r * > 0 with the above property. 

The derivation of this estimate requires not only 
a precise application of an implicit function theorem, 
but also sufficiently accurate knowledge of the general 
solution of equation ( 6) for initial values near , 
rj'j , (j - 1,2) and a time range at least as large as 
the anticipated period T near T * , so that sharp 
estimates of the perturbation of w(t) from Keplerian 
motion w* (t) can be obtained. For this purpose the 
solutions of equation (6) have to be constructed by a 
suitable method of perturbation theory. Because of 
the singularity of equation (6) at w - 0 the motion of 
the satellite is considered in an annulus about the 
planet which contains the precessing Keplerian 
elliptic orbit w* (t) . The value of || P Q || depends 
on the size of this annulus and thus on m, k, e, and 
fx . The main difficulty arises from the fact that 
|| P Q || becomes small along w*(t) only, when a = 
(m/k) V° becomes small, leading to an increase of 
the required range (from 0 to at least 27rk) of the 
independent variable, for which the eccentric or the 
true anomaly can be taken. But thereby the above 
r* decreases with decreasing m/k, almost defeating 
our goal || P Q |[ ^ r* . This difficulty does not appear 
in the classical case of circular w* (t) , or e = 0, 
mentioned earlier. 


IV. CONCLUSION 


Summarizing, we can say that generalization of 
Poincare's small parameter method to the non- 
parametric case by considering the disturbing func- 
tion itself as a generalized small parameter belonging 
to a normed function space leads to applicability of 
classical methods again, and, together with suitable 
periodicity conditions, for example equation (3) , and 
with sufficiently accurate convergent methods of 
perturbation theory, to an existence proof for periodic 
solutions x(t) of equation ( 1) near x* (t) in the 
elliptic satellite case with arbitrary fx in 0 < jx < 1 
[ 2 ]. 
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Finally, by numerical extension of the latter 
families of periodic orbits to greater distances from 
the small planet, interesting new trajectories of the 
restricted three body problem have been found which 
do not belong to the satellite case or to the planetary 
case. Some of these pass repeatedly near both 
attractive bodies [3]; others have been found by my 
collaborator, M. C. Davidson [4], Among the 
latter ones are trajectories which demonstrate the 



phenomenon of temporary capture with satellite 
motion about each one of the attractive bodies and 
periodically alternating transitions from the vicinity 
of one body to the vicinity of the other. Two examples 
( Figs, i and 2) will be included here. They are 
drawn in the rotating x^ x 2 coordinate system under 
the assumption that the attractive bodies represent 
earth E and moon M with \i = 0. 0123 » 1/82. They 
constitute numerical solutions of equation ( 1). 



FIGURE 1 FIGURE 2 
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